Equations are \begin{align} \dot x&=-x^3/3+x-y+I\\ &=-I_{ion}(x,y)+I,\quad I_{ion}(x,y)=x^3-x+y\\ \dot y&=\varepsilon(x-y) \end{align}
using Plots
T=500 #tiempo final, habran 500 tiempos
dt=0.001 #Paso del tiempo
Nt=convert(Int64,T/dt)+1 #Cuantos bines de tiempo son, es decir, se divide el total de tiempos entre el paso dt
t=range(0,T,length=Nt) #Vector de tiempos
function FHNeuler(I,ε,x0,y0) #Función para resolver las ec. FHN con euler
x=zeros(Nt,1)
y=zeros(Nt,1)
x[1]=x0 #Condiciones iniciales
y[1]=y0
for i=2:Nt
x[i]=x[i-1]+dt*(-(x[i-1]^3)/3+x[i-1]-y[i-1]+I) #Ecuacion de x
y[i]=y[i-1]+dt*ε*(x[i-1]-y[i-1]) #Ecuacion en y
end
return x,y
end
I=-0.41 #Corriente aplicada que da una dinámica neuronal estable
ε=0.05
x0=0.0 #randn(1)[1] Condiciones iniciales, pueden ser aleatorias
y0=-1.0 #randn(1)[1]
@time x,y=FHNeuler(I,ε,x0,y0); #llamamos a la función y guardamos el output en variables x y
#Plot de los vectores solucion de las ecuaciones FHN, dado I=-0.41
#plot(t,x,label="x", linewidth=4, ylabel="Voltage (mV)",xlabel="Time (ms)", title= "FitzHugh-Nagumo model")
#plot!(t,y,label="y", linewidth= 2) #El signo de admiración indica que las gráficas estarán juntas
#Comente esta parte del código porque enseguida agrego la figura
using DifferentialEquations
function FHNdiffEq(du,u,p,t)
I=p[1]
ε=p[2]
x=u[1]
y=u[2]
du[1]=-x^3/3+x-y+I
du[2]=ε*(x-y)
end
T=500.0
tspan=(0.0,T)
I=-0.41
ε=0.05
x0=0.0 #randn(1)[1]
y0=-1.0 #randn(1)[1]
p=(I,ε)
prob = ODEProblem(FHNdiffEq,[x0;y0],tspan,p)
@time sol = solve(prob,reltol=10e-9,abstol=10e-9,AutoTsit5(Rosenbrock23())); ## See Julia Diff EQ page
p1= plot(sol, linewidth= 2, xlabel="Time (ms)", ylabel="Voltage (mV)")
p2= plot(sol[1,:],sol[2,:],label="(x(t),y(t))",xlabel="x",ylabel="y", linewidth=2)
plot(p1,p2, layout=(1,2), size=(750,350))
x = -3:0.1:3
y = -3:0.1:3
Xnullcline(x,y)=-x^3/3+x-y+I
Ynullcline(x,y)=x-y
plot(sol,vars=(1,2),label="(x(t),y(t))",xlabel="x",ylabel="y", xlim=(-3.0, 3), linewidth= 3)
contour!(x,y,Xnullcline,levels=[0])
contour!(x,y,Ynullcline,levels=[0],colorbar=false, color= "green")
#FHN phase portraits with different initial values for x and y
x = -3:0.01:3
y = -3:0.01:3
Xnullcline(x,y)=-x^3/3+x-y+I
Ynullcline(x,y)=x-y
contour(x,y,Xnullcline,levels=[0])
contour!(x,y,Ynullcline,levels=[0],colorbar=false)
for x0 in range(-3,3,length=10)
for y0 in range(-3,3,length=10)
prob = ODEProblem(FHNdiffEq,[x0;y0],tspan,p)
sol = solve(prob)
plot!(sol[1,:],sol[2,:],lw=0.5,legend=false)
end
end
plot!()
Approximate the value of $I_{Hopf}$ at which the model undergoes a Hopf bifurcation. How does $I_{Hopf}$ depend on $\varepsilon$? For instance, does $I_{Hopf}$ increase or decrease as $\varepsilon$ increases? Use simulations and, possibly, linear algebra to justify your answer.
A Hopf bifurcation is a critical point where a system's stability switches and a periodic solution arises.
The appearance or the disappearance of a periodic orbit through a local change in the stability properties of a fixed point is known as the Hopf bifurcation.
In the Andronov-Hopf bifurcation the resting state in a neuronal model does not disappear, but it loses stability.
This loss of stability is accompanied by one of the following two dynamics:
1) The apperance of a stable limit cycle appearing from a stable equilibrium (i.e. supercritical Andronov Hopf).
2) The disapperance of an unstable limit cycle to a stable fixed equilibrium (i.e. subcritical Andronov Hopf).
The next bifurcation diagram corresponds to the Supercritical (left) or Subcritical (right) Andronov-Hopf bifurcation, taken from Izhikevich, E. 2007.
A neuron is excitable because its resting state is near a bifurcation, i.e. transition from resting to periodic spiking. Such bifurcation (transition) can be revealed by applying a ramp of current.
The following figure shows the transition from the resting state to repetitive spiking via Andronov-Hopf bifurcation. Top figure is taken from Izhikevich, 2007. Bottom figure shows the dynamics that I found simulating a current ramp in the FHN model.
In this section I simulate the supercritical Hopf bifurcation.
I focused on this type of Hopf bifurcation because according to the literature it is the one that appears most commonly in neurons and their models. Furthermore, subcritical Hopf is reported to be more difficult and unusual to find.
Through various simulations and based on previous literature, I found that $I_{Hopf} = -0.3147$ in the FitzHugh Nagumo model.
T=500.0
tspan=(0.0,T)
#Aquà voy a establecer 5 simulaciones, la I1 corresponde a la bifurcación de HOPF Super, I4 al resting state del modelo, previamente visto.
#Las otras tres simulaciones corresponden a la transicion del resting a la bifurcacion, esto se logra disminuyendo la corriente aplicada.
#En I=0 estamos en el estado oscilatorio del sistema
I0= 0.0 #Estado oscilatorio
I1= -0.3147 #HOPF super
I2= -0.31475
I3= -0.33
I4= -0.41 #Resting State
#Es importante notar que el sistema FHN es SIMÉTRICO, asà que si cambiamos los signos de I1 a I5, observaremos el
#mismo comportamiento solo que en el dominio positivo.
#Me dediqué a simular en el dominio negativo porque es fisiológicamente donde tiene sentido el comportamiento.
#Claro, pensando que el potencial de membrana de la neurona oscila entre -40 y -60 mV.
ε=0.05
x0=0.0
y0=-1.0
p0=(I0,ε)
p1=(I1,ε)
p2=(I2,ε)
p3=(I3,ε)
p4=(I4,ε)
prob0 = ODEProblem(FHNdiffEq,[x0;y0],tspan,p0)
prob1 = ODEProblem(FHNdiffEq,[x0;y0],tspan,p1)
prob2 = ODEProblem(FHNdiffEq,[x0;y0],tspan,p2)
prob3 = ODEProblem(FHNdiffEq,[x0;y0],tspan,p3)
prob4 = ODEProblem(FHNdiffEq,[x0;y0],tspan,p4)
@time sol0 = solve(prob0,reltol=10e-9,abstol=10e-9,AutoTsit5(Rosenbrock23()));
@time sol1 = solve(prob1,reltol=10e-9,abstol=10e-9,AutoTsit5(Rosenbrock23()));
@time sol2 = solve(prob2,reltol=10e-9,abstol=10e-9,AutoTsit5(Rosenbrock23()));
@time sol3 = solve(prob3,reltol=10e-9,abstol=10e-9,AutoTsit5(Rosenbrock23()));
@time sol4 = solve(prob4,reltol=10e-9,abstol=10e-9,AutoTsit5(Rosenbrock23()));
#Iosc=0.0
x = -3:0.1:3
y = -3:0.1:3
Xnullcline(x,y)=-x^3/3+x-y+I0
Ynullcline(x,y)=x-y
p1= plot(sol0, linewidth= 2, ylabel="Voltage (mV)", xlabel= "Time (ms)")
p2=plot(sol0[1,:],sol0[2,:],label="(x(t),y(t))",xlabel="x",ylabel="y")
p3=plot(sol0,vars=(1,2),label="(x(t),y(t))",xlabel="x",ylabel="y", xlim=(-3.0, 3), linewidth=2, size=(500,400))
p3=contour!(x,y,Xnullcline,levels=[0])
p3=contour!(x,y,Ynullcline,levels=[0],colorbar=false, color="green")
plot(p1, p2, p3, layout=(1,3), legend= false,size=(980,300), title="Iosc = 0.0")
#IHopf = -0.3147
x = -3:0.1:3
y = -3:0.1:3
Xnullcline(x,y)=-x^3/3+x-y+I1
Ynullcline(x,y)=x-y
p1= plot(sol1, linewidth= 2, ylabel="Voltage (mV)", xlabel= "Time (ms)")
p2=plot(sol1[1,:],sol1[2,:],label="(x(t),y(t))",xlabel="x",ylabel="y")
p3=plot(sol1,vars=(1,2),label="(x(t),y(t))",xlabel="x",ylabel="y", xlim=(-3.0, 3), linewidth=2, size=(500,400))
p3=contour!(x,y,Xnullcline,levels=[0])
p3=contour!(x,y,Ynullcline,levels=[0],colorbar=false, color="green")
plot(p1, p2, p3, layout=(1,3), legend= false,size=(980,300), title= "IHopf = -0.3147")
#I = -0.31475
x = -3:0.1:3
y = -3:0.1:3
Xnullcline(x,y)=-x^3/3+x-y+I2
Ynullcline(x,y)=x-y
p1= plot(sol2, linewidth= 2, ylabel="Voltage (mV)", xlabel= "Time (ms)")
p2=plot(sol2[1,:],sol2[2,:],label="(x(t),y(t))",xlabel="x",ylabel="y")
p3=plot(sol2,vars=(1,2),label="(x(t),y(t))",xlabel="x",ylabel="y", xlim=(-3.0, 3), linewidth=2, size=(500,400))
p3=contour!(x,y,Xnullcline,levels=[0])
p3=contour!(x,y,Ynullcline,levels=[0],colorbar=false, color="green")
plot(p1, p2, p3, layout=(1,3), legend= false, size=(980,300), title= "I = -0.31475")
#I= -0.33
x = -3:0.1:3
y = -3:0.1:3
Xnullcline(x,y)=-x^3/3+x-y+I3
Ynullcline(x,y)=x-y
p1= plot(sol3, linewidth= 2, ylabel="Voltage (mV)", xlabel= "Time (ms)")
p2=plot(sol3[1,:],sol3[2,:],label="(x(t),y(t))",xlabel="x",ylabel="y")
p3= plot(sol3,vars=(1,2),label="(x(t),y(t))",xlabel="x",ylabel="y", xlim=(-3.0, 3), linewidth=2, size=(500,400))
p3=contour!(x,y,Xnullcline,levels=[0])
p3=contour!(x,y,Ynullcline,levels=[0],colorbar=false, color="green")
plot(p1, p2,p3, layout=(1,3), legend= false,size=(980,300), title= "I = -0.33")
#I RS= -0.41
x = -3:0.1:3
y = -3:0.1:3
Xnullcline(x,y)=-x^3/3+x-y+I4
Ynullcline(x,y)=x-y
p1= plot(sol4, linewidth= 2, ylabel="Voltage (mV)", xlabel= "Time (ms)")
p2=plot(sol4[1,:],sol4[2,:],label="(x(t),y(t))",xlabel="x",ylabel="y")
p3=plot(sol4,vars=(1,2),label="(x(t),y(t))",xlabel="x",ylabel="y", xlim=(-3.0, 3), linewidth=2, size=(500,400))
p3=contour!(x,y,Xnullcline,levels=[0])
p3=contour!(x,y,Ynullcline,levels=[0],colorbar=false, color="green")
plot(p1, p2,p3, layout=(1,3), legend= false, size=(980,300), title= "IRS = -0.41")
The following figure shows the transition of the FitzHugh Nagumo model from the Resting State ($I_{RS} = -0.6$) to the oscillatory state ($I_{OSC} = 0.0$) passing through the Supercritical Hopf bifurcation ($I_{Hopf} = -0.3147$). As the bifurcation parameter $I$ increases, the stable equilibrium loses stability and gives birth to a stable limit cycle with growing amplitude. Also shown is the bifurcation diagram of the simulated process (bottom right), taken from Izhikevich, E. 2007.
T=500.0
tspan=(0.0,T)
I=-0.3147 #Hopf bifurcation I parameter
x0= 0.0
y0= -1.0
x = -3:0.01:3
y = -3:0.01:3
Xnullcline(x,y)=-x^3/3+x-y+I
Ynullcline(x,y)=x-y
contour(x,y,Xnullcline,levels=[0])
contour!(x,y,Ynullcline,levels=[0],colorbar=false, color="green")
for ε in range(0.001,1,length=10)
p= (I,ε)
prob = ODEProblem(FHNdiffEq,[x0;y0],tspan,p)
sol = solve(prob)
plot!(sol[1,:],sol[2,:],lw=2, label= ε, title="FHN model I Hopf with variable epsilon", xlabel="x", ylabel="y")
end
plot!()
Dynamics of the FHN model at the Hopf bifurcation with different values of $\varepsilon$
First, we must remember that $\varepsilon$ is found in the FHN model in the following equation $\dot y=\varepsilon(x-y)$. This tells us that the value of epsilon acts as a "scaler" of the equation $\dot y$, that is, its value causes the equation to increase or decrease, which translates into fast or slow dynamics of the variable $y$ in the model, respectively.
Therefore, the value of $\varepsilon$ is important because it dictates the kinetics of the recovery variable $y$, and therefore the dynamics of the system.
With this in mind, we can now analyze what happens to the system dynamics at the Hopf bifurcation, and thus find how the value $I_{Hopf}$ depends on $\varepsilon$.
#3 more simulations, one with a very small epsilon (0.001), a "normal" epsilon (0.05) and lastly, with a big epislon value (5)
T=500.0
tspan=(0.0,T)
I=-0.3147 #Hopf bifurcation I parameter
x0= 0.0
y0= -1.0
ε1= 0.001
ε2= 0.05
ε3= 5
p1= (I,ε1)
p2= (I,ε2)
p3= (I,ε3)
prob1 = ODEProblem(FHNdiffEq,[x0;y0],tspan,p1)
prob2 = ODEProblem(FHNdiffEq,[x0;y0],tspan,p2)
prob3 = ODEProblem(FHNdiffEq,[x0;y0],tspan,p3)
sol1 = solve(prob1)
sol2 = solve(prob2)
sol3 = solve(prob3)
#Plot phase portrait
x = -3:0.01:3
y = -3:0.01:3
Xnullcline(x,y)=-x^3/3+x-y+I
Ynullcline(x,y)=x-y
plot(sol1[1,:],sol1[2,:],lw=3, label= ε1, title="FHN model I Hopf with 3 values of epsilon", xlabel="x", ylabel="y", color= "#ffa41b")
plot!(sol2[1,:],sol2[2,:], lw=1, label= ε2, color="#00bdaa", ls= :dash)
plot!(sol3[1,:],sol3[2,:],lw=2, label= ε3, color="#fe346e")
contour!(x,y,Xnullcline,levels=[0], color="black")
contour!(x,y,Ynullcline,levels=[0],colorbar=false, color="black")
Dynamics of the FHN model at the Hopf bifurcation given 3 values of $\varepsilon$.
In this figure we can easily observe the dependence of $I_{Hopf}$ with more "notable" $\varepsilon$ values.
When $\varepsilon = 0.001$ (yellow), the system is truncated at the zerocline of $x$, as in the previous graph, and therefore, the Hopf bifurcation does not exist.
When $\varepsilon = 0.05$ (blue), the Hopf bifurcation occurs at $I = -0.3147$, as we have been seeing throughout the exercise.
Finally, when $\varepsilon = 5$ (pink), the bifurcation disappears and the trajectory of the system from the initial conditions to the stable fixed point is almost instantaneous. If we wanted to retrieve $I_{Hopf}$, the applied current $I >>> -0.3147$.
Modify the code to be able to apply steps of applied currents and train of current pulses with variable amplitude, frequency and duty cycle. Simulate the model for current steps of various amplitude and for pulse trains of various amplitude, frequency, and duty cycle. Describe the observed behavior.
#Función para el pulso (step) de corriente
STEP(t,tstep,Istep)=(sign(t-tstep)+1)/2*Istep
T=500
dt=0.001
Nt=convert(Int64,T/dt)+1
t=range(0,T,length=Nt)
function FHNeulerSTEP(I, ε, x0, y0, tstep, Istep)
"""Función para resolver las ec. FHN con euler,
con la capacidad de aplicar pulsos de corriente (steps)
Argumentos:
I= corriente aplicada inicial
ε= variable que controla la "velocidad" de evolución de y
x0= condición inicial del sistema en x
y0= condición inicial del sistema en y
tstep= tiempo en el cual se dará el pulso de corriente
Istep= corriente del pulso
"""
x=zeros(Nt,1)
y=zeros(Nt,1)
x[1]=x0
y[1]=y0
for i=2:Nt
tin= i*dt
x[i]=x[i-1]+dt*(-(x[i-1]^3)/3+x[i-1]-y[i-1]+I+ STEP(tin, tstep,Istep))
y[i]=y[i-1]+dt*ε*(x[i-1]-y[i-1])
end
return x,y
end
T=500.0
tspan=(0.0,T)
tstep= 250.0
Istep= -1.0 #Corriente del pulso
I= -0.45 #Aquà el sistema esta en un "edo estable", para la Bifurcación de Hopf +- 0.3147319
ε= 0.05
x0= 0.0
y0= -1.0
p=(I,ε)
@time x,y= FHNeulerSTEP(I, ε, x0, y0, tstep, Istep);
#plot1= plot(t,x,label="x", linewidth=3, title="I step= -1.0")
#plot1= plot!(t,y,label="y", linewidth=3)
#plot2= plot(t, I.+STEP.(t,tstep,Istep), linewidth=3, ylabel="I app", xlabel="Time (ms)")
#plot(plot1, plot2, layout=(2,1))
#Comenté esta parte del código porque enseguida agrego las gráficas obtenidas mediante diferentes simulaciones
Simulation of current steps applied to the neuronal model of FHN. The initial condition of the system is $I = -0.45$ in all cases. The step arrival time is $t_{step} = 250 ms$; the only parameter that changes is the applied step current ($I_{step}$).
Note: The y axis of the "spike" graph represents Voltage (mV).
We can see that in the case of $I_{step} = -0.45$ (upper left), at the beginning we see a spike due to the initial conditions of the system (as in all cases), and then a hyperpolarization upon arrival of the pulse, that is, a decrease in the membrane voltage. This is because we went from being in $I_{app} = -0.45$ to $I_{app} = -0.9$.
In the graph of $I_{step} = 0.45$ (upper right), after the resting state and at the arrival of the step, sustained spikes are generated throughout the pulse. This occurs by adding positive current to the system that pushes it to depolarization.
In the graph of $I_{step} = -1.0$ (lower left), we observe the same effect as in $I_{step} = -0.45$, but the hyperpolarization is slightly higher, this because we bring the membrane potential (mV) to more negative values. This happens because we go from $I_{app} = -0.45$ (initial conditions) to $I_{app} = -1.4$.
Finally, in the graph of $I_{step} = 1.0$ (lower right), we can observe that when the step arrives there is an initial depolarization that, because the injected current is large and sustained, is maintained over time. In this case, the system cannot be repolarized and loses its resting state.
###Función para el tren de pulsos
function TRAIN(t, freq, Istep, DC)
if DC==0
train=(sign(sin(2*pi*freq*t)*DC)/2)*Istep #Si el DC =0, que todo se multiplique por 0, ya que significa que no existen pulsos
elseif DC==1
train= (sign(t)+1)/2*Istep #Cuando el DC es 1, significa que el step sube en el inicio y se queda arriba todo el experimento
elseif DC<0.5
train= (sign(sin(2*pi*freq*t)+(DC-1))/2)*Istep #Cuando el valor DC<0.5, hay que restarle 1 porque p.e. si ponemos DC=0.2 el tren lo interepreta como DC=0.8
elseif DC>0.5
train= (sign(sin(2*pi*freq*t)+(DC))/2)*Istep #Cuando el valor del DC >0.5, se queda igual
elseif DC==0.5
train=(sign(sin(2*pi*freq*t))/2)*Istep #Si DC=0.5, entonces no se le debe agregar nada al seno, ya que el tren sera periodico, ie estara 0.5 abajo y 0.5 de tiempo arriba
end
return train
end
T=500
dt=0.001
Nt=convert(Int64,T/dt)+1
t=range(0,T,length=Nt)
function FHNeulerTRAIN(I, ε, x0, y0, freq, Istep, DC)
"""Función para resolver las ec. FHN con euler,
con la capacidad de aplicar trenes de pulsos de corriente
Argumentos:
I= corriente aplicada inicial
ε= variable que controla la "velocidad" de evolución de y
x0= condición inicial del sistema en x
y0= condición inicial del sistema en y
tstep= tiempo en el cual se dará el pulso de corriente
Istep= corriente del pulso
"""
x=zeros(Nt,1)
y=zeros(Nt,1)
x[1]=x0
y[1]=y0
for i=2:Nt
tin= i*dt
x[i]=x[i-1]+dt*(-(x[i-1]^3)/3+x[i-1]-y[i-1]+ I+TRAIN(tin, freq, Istep, DC))
y[i]=y[i-1]+dt*ε*(x[i-1]-y[i-1])
end
return x,y
end
T=500.0
tspan=(0.0,T)
Istep= 0.3 #Corriente del pulso
I= -0.45 #En 0.3147 esta la bifurcacion de Hopf
ε= 0.05
x0= 0.0
y0= -1.0
freq= 0.01
DC= 0.5
p=(I,ε)
@time x,y= FHNeulerTRAIN(I, ε, x0, y0, freq, Istep, DC);
#plot1= plot(t,x,label="x", linewidth=3)
#plot1= plot!(t,y,label="y", linewidth=3)
#plot2= plot(t, I.+TRAIN.(t, freq, Istep, DC), linewidth=3, title="DC=0.5, freq=0.01, I step=0.3", xlabel="Time(ms)", ylabel="I app")
#plot(plot1, plot2, layout=(2,1))
#Comenté esta parte del código porque enseguida agrego las gráficas obtenidas mediante diferentes simulaciones
Simulations of current pulse trains applied to the FHN neural model. The base current is $I = -0.45$ in all cases .
Use phase portrait analysis to explain why the model responds with a rebound spike to negative current pulses.
The question that gave rise to this phenomenon is that if excitatory inputs depolarize the membrane potential, ie bring the neuron closer to its firing threshold, and the inhibitory inputs hyperpolarize the potential and move it away from said threshold, then how is it neurons can fire in response to inhibitory inputs?
Example of the described phenomena in a rat's brainstem neuron:
This phenomenon is observd in th FHN model and it is called anodal break excitation, rebound spike or postinhibitory spike.
Many biologists say that rebound responses are due to the activation and inactivation of certain slow currents, which bring the membrane potential over the threshold. The problem with this explanation is that neither the FHN model nor the neuron in the above figure has these currents, and even if they did, the hyperpolarization is too short and too weak to affect them.
The following GIF shows the Anodal break excitation (post-inhibitory rebound spike) in the FitzHugh-Nagumo model, taken from Eugene M. Izhikevich and Richard FitzHugh (2006), Scholarpedia, 1(9):1349.
As the stimulus $I$ becomes negative (hyperpolarization), the resting state shifts to the left. As the system is released from hyperpolarization, the trajectory starts from a point far below the resting state, makes a large-amplitude excursion, i.e. fires a transient spike, and then returns to the resting state.
T=500.0
tspan=(0.0,T)
Istep= 0.3 #Corriente del pulso
I= -0.45 #En 0.3147 esta la bifurcacion de Hopf
ε= 0.05
x0= 0.0
y0= -1.0
freq= 0.01
DC= 0.5
p=(I,ε)
@time x,y= FHNeulerTRAIN(I, ε, x0, y0, freq, Istep, DC);
X = -3:0.1:3
Y = -3:0.1:3
Xnullcline(X,Y)=-X^3/3+X-Y+I
Xnullcline2(X,Y)=-X^3/3+X-Y+(I-Istep)
Ynullcline(X,Y)=X-Y
p1= plot(t,x, linewidth=3, xlabel="Time (ms)", ylabel="Voltage (mV)", label=false, xlim= (228,340))
p2=plot(x, y, lw=2, xlabel="x", ylabel="y", label=false)
p3=plot(x, y,label="(x(t),y(t))",xlabel="x",ylabel="y", xlim=(-3.0, 3), linewidth= 3)
p3=contour!(X,Y,Xnullcline,levels=[0])
p3=contour!(X,Y,Xnullcline2,levels=[0], ls= :dash)
p3=contour!(X,Y,Ynullcline,levels=[0],colorbar=false, color= "green")
plot(p1,p2,p3, layout=(1,3), size=(980,320))
Simulation of the rebound spike phenomenon in the FHN model.
The first figure (left) is the simulation of the rebound spike in response to a short hyperpolarizing pulse. The figure in the center corresponds to the evolution of the variables $x$ and $y$ over time. Finally, the third figure shows the phase portrait of the FHN model: the dynamic $x(t), y(t)$ is shown in blue, the solid red nullcline corresponds to $I = -0.45$, and the dotted red nullcline corresponds to $I = -0.6$.
Phase portrait analysis.
Our system starts at the solid nullcline (original resting state, $I = -0.45$), then the hyperpolarizing stimulus arrives (step up), causing the resting state to shift to the left, i.e. brings the system to the dashed nullcline ($I = -0.6$). When the hyperpolarizing pulse is removed (step down), the trajectory starts in the new resting state ($-0.6$), then travels the orbit corresponding to a spike and returns to its original resting state ($-0.45$). In this case, since the simulation was done with a train of pulses, several trajectories are seen in the phase portrait.
recta(pendiente, t, tstart)=pendiente*(t - tstart)
#recta es la función que permite modificar la velocidad de la rampa, através del argumento "pendiente"
RAMP(t, tstart, pendiente)=STEP(t, tstart, 1.0)* recta(pendiente, t, tstart)
T= 1300.0
dt=0.001
Nt=convert(Int64,T/dt)+1
t=range(0,T,length=Nt)
function FHNeulerRAMP(I, ε, x0, y0, tstart, pendiente)
"""Función para resolver las ec. FHN con euler,
con la capacidad de aplicar rampas de corriente
Argumentos:
I= corriente aplicada inicial
ε= variable que controla la "velocidad" de evolución de y
x0= condición inicial del sistema en x
y0= condición inicial del sistema en y
tstart= el tiempo en el que se da la rampa, es equivalente a tstep
pendiente= valor de la pendiente de la rampa
"""
x=zeros(Nt,1)
y=zeros(Nt,1)
x[1]=x0
y[1]=y0
for i=2:Nt
tin= i*dt
x[i]=x[i-1]+dt*(-(x[i-1]^3)/3+x[i-1]-y[i-1]+I+ RAMP(tin, tstart, pendiente))
y[i]=y[i-1]+dt*ε*(x[i-1]-y[i-1])
end
return x,y
end
T=1300.0
tstart= 50.0
pendiente= 0.001
I= -0.45
ε= 0.05
x0= 0.0
y0= -1.0
p=(I,ε)
@time x,y= FHNeulerRAMP(I, ε, x0, y0, tstart, pendiente);
#plot1= plot(t,x,label="x", linewidth=3)
#plot1= plot!(t,y,label="y", linewidth=3)
#plot2= plot(t, I.+RAMP.(t, tstart, pendiente), linewidth=3, title="m= 0.001, I start= -0.6", xlabel="Time (ms)", ylabel= "I app" )
#plot(plot1, plot2, layout=(2,1))
#Comenté esta parte del código porque enseguida agrego las gráficas obtenidas mediante diferentes simulaciones
The excitation block is a phenomenon presented in many neurons in a physiological way. It consists of the cessation of repetitive spiking as the amplitude of the applied current increases.
The excitation block phenomenon in a rat's neuron and in a simulation of the FHN model. The top figure shows that the spiking activity of the rat's neuron is blocked by strong excitation (i.e., injecting a strong depolarizing current), taken form Izhikevich, 2007. The figure below shows the same phenomenon in my simulation of the FHN model with an ramp of injected current. Ramp parameters: $m=0.001$, $I_{start} = -0.6$.
X = -3:0.1:3
Y = -3:0.1:3
#Xnullcline(X,Y)=-X^3/3+X-Y+I Esta es la que empieza en -0.6, pero preferà quitarla por estética
Xnullcline2(X,Y)=-X^3/3+X-Y+0
Xnullcline3(X,Y)=-X^3/3+X-Y+0.3147
Xnullcline4(X,Y)=-X^3/3+X-Y-0.3147
Ynullcline(X,Y)=X-Y
plot(x,y, linewidth=2, xlabel="x", ylabel="y", label=false)
contour!(X,Y,Xnullcline,levels=[0])
contour!(X,Y,Xnullcline2,levels=[0], color="orange")
contour!(X,Y,Xnullcline3,levels=[0], color="black")
contour!(X,Y,Xnullcline4,levels=[0], color="purple")
contour!(X,Y,Ynullcline,levels=[0],colorbar=false, color="green")
Phase portrait analysis of the excitation block phenomenon in the FHN model.
The figure on the left shows the phase portrait corresponding to the simulation of the excitation block of the FHN model with a ramp of applied current shown in the previous cell. The gif on the right shows the phase portrait of the excitation block through time, taken from Eugene M. Izhikevich and Richard FitzHugh (2006), Scholarpedia, 1(9):1349.
| |
When the current $I$ is low ($I_{RS} = -0.45$), the equilibrium (intersection of nullclines) is in the left arm of the x-nullcline (stable), therefore the system is in its resting state. Increasing $I$ moves the x-nullcline upwards reaching the supercritical Hopf bifurcation ($I_{Hopf} = -0.3147$, purple), i.e., the appearance of a stable limit cycle from the disappearance of a stable equilibrium.
If we continue to increase the current ($I_{Osc} = 0.0$, orange) and therefore keep pushing the nullcline upward, the equilibrium slides towards the middle branch (unstable). In this case the model exhibits oscillations, ie, periodic spiking. Increasing the applied current further moves the equilibrium to the right branch of the nullcline (stable), and the spiking activity is blocked (by excitation). The disappearance of this limit cycle occurs via another supercritical Hopf bifurcation ($I_{Hopf} = 0.3147$, black).
In this way, the FHN model can exhibit two Hopf bifurcations in response to ramping up of the injected current, the one leading to the apperance of periodic spiking activity ($I_{Hopf} = -0.3147$), and the one leading to its disappearance ($I_{Hopf} = 0.3147$).
(Optional) Modify the code to run stochastic simulations. Suggestion: Modify the Euler version.
To add "white noise" in a neural model, the Gaussian random number function is used, that is, with distribution centered at 0 and variance of 1.
function FHNeulerStochastic(Dx,Dy,I,ε,x0,y0)
x=zeros(Nt,1)
y=zeros(Nt,1)
x[1]=x0
y[1]=y0
for i=2:Nt
x[i]=x[i-1]+dt*(-(x[i-1]^3)/3+x[i-1]-y[i-1]+I+Dx*randn(1)[1])
y[i]=y[i-1]+dt*ε*(x[i-1]-y[i-1]+Dy*randn(1)[1])
end
return x,y
end
I=-0.45
ε=0.05
x0=0.0 #randn(1)[1] Condiciones iniciales random
y0=-1.0 #randn(1)[1]
Dx= 4.0 #Con Dx > 1.8 en FHN hay umbral y se generan spikes de manera aleatoria
Dy=0.5
@time x,y=FHNeulerStochastic(Dx,Dy,I,ε,x0,y0);
#plot(t,x,label="x", xlabel="Time (ms)", title="Stochastic simulation, Dx= 1.2 and Dy= 0.5")
#plot!(t,y,label="y")
#Comenté esta parte del código porque enseguida agrego las gráficas obtenidas mediante diferentes simulaciones